%%  
% --------------------------------------------------------------
%   1-D, 1-G DIFFUSION-BASED, NONLINEAR RESPONSE MATRIX SOLVER
% --------------------------------------------------------------
%  This test case aims to determine keff for a one-dimensional
%  cartesian slab reactor using analytic response functions that
%  are essentially combinations of transmission/reflection 
%  coefficients and multiples of integrated fluxes.  The current
%  coupling depends implicitly on keff, and a nonlinear problem
%  arises, for which a solution is sought using Newton's method.
%  Additionally, the standard red-black sweeping method has been
%  added for comparison of convergence.  (Hypothesis: nonlinear
%  should be roughly quadratic convergence; sweeping, probably
%  more like linear).
%  j. roberts, 11/2009
clear; clc;

% declare global problem data
global ne dc sa ns sw bcL bcR

problem = 3;

if problem == 1
%------------------------------------------------------------------
% SAMPLE PROBLEM  # 1                                
%  - two slabs of the same material 
%  - zero re-entrant current
% FDM solution yields k ~ 1.002; here, we get 1.0022
ne = 2;                     % number of elements
dc = [ 4.000  4.000 ];      % diffusion coefficient      [cm^2]
sa = [ 0.100  0.100 ];      % absorption cross-section   [1/cm^2]
ns = [ 0.200  0.200 ];      % nu * fission cross-section [1/cm^2]
sw = [ 4.247  4.247 ];      % slab width                 [cm]
bcL = 1;                    % left boundary condition
bcR = 1;                    % right boundary condition
k = 1.0;                    % initial guess of keff
j = ones(1,ne*2);           % initial guess for currents
%------------------------------------------------------------------
elseif problem == 2
%------------------------------------------------------------------
% SAMPLE PROBLEM  # 2                                
%  - five element, infinite core, |mod|f1|mod|f2|mod|
%  - refected bc's
% FDM solution yields k ~ 1.0104; here, we get k ~ 1.0102
ne = 5;
dc = [ 0.200  0.650  0.200  0.600  0.200  ];      
sa = [ 0.020  0.140  0.020  0.150  0.020  ];      
ns = [ 0.000  0.175  0.000  0.170  0.000  ];      
sw = [ 2.000  2.500  4.000  2.500  2.000  ];  
bcL = 1;                    
bcR = 0;  
k = 0.50;
j = ones(1,ne*2);
%------------------------------------------------------------------
elseif problem == 3
% SAMPLE PROBLEM # 3
%  - twenty element, finite core, something asymmetrical 
ne = 80;
dc = [ 0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
       0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
       0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ... 
       0.240  0.660  0.210 1.500  0.510  0.650  0.620  0.630  0.630  0.200  0.200  0.640  0.240  0.600  0.610  0.200  0.620  0.220  0.600  0.200 ];
sa = [ 0.020  0.140  0.020  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
       0.020  0.140  0.020  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
       0.020  0.140  0.020  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020  ...
       0.020  0.140  0.020  0.150  0.120  0.120  0.140  0.120  0.160  0.020  0.020  0.140  0.020  0.100  0.020  0.020  0.140  0.020  0.150  0.020 ];
ns = [ 0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...  
       0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...
       0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ...
       0.000  0.075  0.000  0.090  0.000  0.000  0.081  0.100  0.130  0.000  0.003  0.175  0.060  0.100  0.140  0.000  0.175  0.000  0.370  0.000 ];  
sw = [ 2.200  2.500  4.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ...  
       2.200  2.500  4.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ... 
       2.200  2.500  4.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ... 
       2.200  2.500  4.000  2.500  2.000  2.200  2.500  4.000  2.500  5.000  3.000  2.500  4.000  2.500  2.000  2.000  2.500  1.000  2.500  5.000 ];
bcL = 1;                    
bcR = 1;  
k = 1.9;
j = ones(1,ne*2);   
end

%------------------------------------------------------------------
%---------------------- BEGIN SOLVERS -----------------------------
%------------------------------------------------------------------

% normalize initial currents and form vector of unknowns
j       = j/sum(j);
x       = [j k]';
%xx(1,:) = x;

% -- code for using more robust newton solver
% tol         = [1.d-6, 1.d-6]; 
% params      = [40, 1, 0];
% [result, errs, it_hist] = nsol(x, 'respfct', tol, params);


%---------------------- RED BLACK SOLVER --------------------------
% this solver is not set up optimally from a memory standpoint,
% meaning there are lots of zeros that could be eliminated if
% i wanted to be clever.
tic
[jj,nkrb,itrb] = redblack(ne, dc, sa, ns, sw, x);
t1 = toc;

%---------------------- NONLINEAR SOLVER --------------------------
% the nonlinear function is   
%   -------------------------------     -----
%  | R(k_n)*M - I  |       0       |   | j_n |   
%  |    F(k_n)     | -L(k_n,j_n-1) | x | k_n | = 0
%   -------------------------------     -----
%  where F is 'gains' and L is 'losses'
tic
z       = respfct(x); % initial residual
ro      = norm(z); 
itmx    = 40; 
it      = 0;
zz(1)   = norm(z);
k1      = k; k0 = 0;

%% GMRES %%
tol = [1e-12;1e-12];
[x, it_hist, ierr] = nsolgm(x,'respfct',tol);
t2 = toc;
%%%%%%%%%%%
disp('REDBLACK')
krb = jj(end)
disp('NONLINEAR')
x(1:end-1) = x(1:end-1)/sum(x(1:end-1));
knl = x(end)
disp('COMP')
percdiff = 100*(jj(end)-x(end))/jj(end)
disp('TIMES')
t1,t2
%%
return

while (abs(k1-k0) > 1e-10)% && norm(z) > 1d-10*ro + 1d-10 && it < itmx
    % compute F'(x)
    %fp          = jacob('respfct',x,z)
    %fp          = diffjac(x,'respfct',z)
    % solve F'*s = -F(x); note fp'fp for spd
    s           = -(fp'*fp)\(fp'*z);
    %[q,r] = qr(fp);
    %s = r\(q'*z);
    % create new x
    k0          = x(end);
    x           = x + s;
    k1          = x(end);
    % normalize currents
    x(1:end-1)  = x(1:end-1)/sum(x(1:end-1));
    it          = it+1;
%     xx(it,:)  = x;
    z           = respfct(x);
    zz(it)      = norm(z);
    nknl(it)    = norm(k1-k0);   
end
t2 = toc;

%--- OUTPUT
 
semilogy(1:itrb,nkrb,'k',1:it,nknl,'r--',1:.1:max(itrb,it),1e-10,'k--')
xlabel('iteration')
ylabel('residual in k')
legend('red-black','newton')

  disp(' ')
  disp([' *** final results ***'])
  disp(' ')
  disp(['       |    red-black    |    nonlinear    |   '])
fprintf('    it |      %3i        |      %3i        | \n', itrb,it)
fprintf('  keff |%16.13f |%16.13f | \n', jj(end),x(end))
p1 = log( nkrb(end) / nkrb(end-1) )/log( nkrb(end-1)/nkrb(end-2));
p2 = log( nknl(end) / nknl(end-1) )/log( nknl(end-1)/nknl(end-2));
fprintf('     p |      %4.3f      |      %4.3f      | \n',p1,p2)
fprintf('  time |      %4.3f      |      %4.3f      | \n',t1,t2)

% this is the jacobian for the first problem with some analytically
% derived values; the matrix is singular to working precision
% AA = [ -1       0       0.7999247629           0       -0.3917587446e-1*x(2) 
%         0       -1      0.2197291477    0       -0.4031917547e-1*x(2)
%         0       -1      0.2197291477    0       -0.4031917547e-1*x(2)
%         0      0.7999247629  0              -1       -0.3917587446e-1*x(2)
%        -0.1824010047   0.8395989952  0.8395989952  -0.182  -2*(1.894426949*x(1)+0.8944269490*x(2))
%        ]
